function [Gamma_New, F_new] = IPCA_base(Gamma_Old,X,W,Nt,varargin)

% Gamma_Old is L*Ktilde matrix
% Ktilde is K+M
% M is the number of pre-specified factors
% If there is only alpha, then M=1 and PSF(t) = 1 for all t

% Gamma_New is L*Ktilde
% F_New is KxT
% If K=0, then F_New = [] and Gamma_New = Gamma_delta

% Standard IPCA call as IPCA_base(Gamma_Old,X,W,Nt) with Gamma_Old as L*K
% IPCA with PSF call as IPCA_base(Gamma_Old,X,W,Nt,PSF) with Gamma_Old as L*(K+M)
% IPCA with alpha call as IPCA_base(Gamma_Old,X,W,Nt,ones(1,T)) with Gamma_Old as L*(K+1)

if nargin>4 && ~isempty(varargin{1})
    PSF = varargin{1};
    PSF_version = true;
    M = size(PSF,1);
    [L,Ktilde] = size(Gamma_Old);
    K = Ktilde - M;
else
    PSF_version = false;
    [L,K] = size(Gamma_Old);
    Ktilde = K;
end

T = length(Nt);

% get F_new
F_new = [];
if K>0
    F_new = NaN(K,T);
    if PSF_version
        for t=1:T-1
            F_new(:,t+1) = (Gamma_Old(:,1:K)'*W(:,:,t)*Gamma_Old(:,1:K) )\( Gamma_Old(:,1:K)'*(X(:,t+1) - W(:,:,t)*Gamma_Old(:,K+1:Ktilde)*PSF(:,t+1)) );
        end
    else
        for t=1:T-1
            F_new(:,t+1) = (Gamma_Old'*W(:,:,t)*Gamma_Old )\( Gamma_Old'*X(:,t+1) );
        end
    end
end

% get Gamma_new
Numer = zeros(L*Ktilde,1);
Denom = zeros(L*Ktilde,L*Ktilde);
if PSF_version
    if K>0
        Ftilde_new = [F_new; PSF];
    else
        Ftilde_new = PSF;
    end
else
    Ftilde_new = F_new;
end
for t=1:T-1
    Numer = Numer + kron( X(:,t+1) , Ftilde_new(:,t+1)                    )*Nt(t+1);
    Denom = Denom + kron( W(:,:,t) , Ftilde_new(:,t+1)*Ftilde_new(:,t+1)' )*Nt(t+1);
end
Gamma_New_trans_vec = Denom\Numer;
Gamma_New_trans     = reshape(Gamma_New_trans_vec,Ktilde,L);
Gamma_New           = Gamma_New_trans';

if K>0
    % GammaBeta orthonormal and F_New Orthogonal
    R1                  = chol(Gamma_New(:,1:K)'*Gamma_New(:,1:K),'upper');
    [R2,~,~]            = svd(R1*F_new(:,2:T)*F_new(:,2:T)'*R1');
    Gamma_New(:,1:K)    = (Gamma_New(:,1:K)/R1)*R2;
    F_new(:,2:T)        = R2\(R1*F_new(:,2:T));

    % Sign convention on GammaBeta and F_New
    sg = sign(mean(F_new(:,2:T),2));
    sg(sg==0)=1; % if mean zero, do not flip signs of anything -- cannot impose this sign convention
    Gamma_New(:,1:K) = Gamma_New(:,1:K).*sg';
    F_new(:,2:T) = F_new(:,2:T) .* sg;
end

% Orthogonality between GammaBeta and GammaDelta
if PSF_version && K>0
    Gammabeta = Gamma_New(:,1:K);
    Gammadelta = Gamma_New(:,K+1:end);
    Gammadelta = ( eye(L) - Gammabeta*Gammabeta' )*Gammadelta;

    gamma = Gammabeta'*Gammadelta; %K*M reg coef
    F_new = F_new + gamma*PSF;

    Gamma_New = [Gammabeta Gammadelta];
    
    % AGAIN: Sign convention on GammaBeta and F_New
    sg = sign(mean(F_new(:,2:T),2));
    sg(sg==0)=1; % if mean zero, do not flip signs of anything
    Gamma_New(:,1:K) = Gamma_New(:,1:K).*sg';
    F_new = F_new .* sg;
end
return